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Abstract 

We propose a statistical formulation to estimate possible errors in long baseline neutrino oscillation 
experiments caused by uncertainties in the Earth matter density. A quantitative investigation of the 
effect is made on the CP asymmetry in future neutrino factory experiments. 

Leptonic CP violation (CPV) is one of the main challenges in future long baseline (LBL) neutrino 
oscillation experiments fl, where more accurate measurements of the neutrino oscillation parameters are 
anticipated. However, since the neutrino beam travels a long path through the Earth, the MSW matter effect 
can mimic a non- vanishing CP phase and makes it non-trivial to extract the intrinsic CP phase. Therefore, 
a thorough quantitative understanding of the matter effect and its possible uncertainties is necessary before 
an accurate account of the CPV effect can be achieved. 

The critical quantity in discerning the matter effect is the value of the Earth matter density as a function 
of the baseline. A number of approaches have been suggested on using different Earth density profiles of the 
Earth for the electron number density (DEN). In these works, DEN is taken either as a distance-averaging 
effective constant [S , an adiabatic approximation profile , mantle-core-mantle layers approximation j^] , 
multi-step functions or the preliminary reference Earth model (PREM) 0. We refer to Ref. || for a 
review of some of the available Earth density models. It is clear from these works that the analysis of the 
leptonic CPV will depend on the Earth density model adopted to analyze the experimental data. 

From the geophysics point of view, a density model is understood to be an approximation of the Earth 
and has its inherent uncertainties. For a brief discussion on the uncertainty of DEN, we refer to Refs. |J and 
||. Detailed discussions can be found in geophysics revie w ffo| , [Tl| , [i"2| . As to PREM, significant deviations 
from PREM due to local variation have been documented Hand averaged per spherical shell with thickness 
of 100 Km, its precision is roughly 5%. 

Instead of examining quantitatively all the schemes of the Earth density to find out the preferred density 
scheme to use, we would like to raise and attempt to provide at least a partial answer to the following basic 
questions: What is the effect of uncertainties in the Earth matter density on the measurement of CPV? Can 
the effect be so severe that it demands a more accurate knowledge of the Earth density before a reliable CP 
phase can be extracted in LBL experiments? 

We propose in this paper to estimate the error due to the uncertainty in the Earth matter density as 
the variance of a CP odd oscillation probability. As described below the procedure in making the estimate 
is a weighted average over the whole sample space of possible Earth density profiles. We first introduce 
a local variance function a(x) to characterize the density uncertainties at each point of the Earth along a 
given baseline. With the variance function we sample the earth density profiles, then construct an averaging 
process to calculate the deviation of a physical quantity from its mean value. The probability of the various 
DEN examples is taken as a logarithmic distribution suitable for non-negative quantities, although other 
statistics approach may be used. 
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We begin the formulation with the flavor Hamiltonian that governs the neutrino propagation in matter. 
Omitting terms that leads to only a common phase to all flavor states, we have, in the scheme of three flavors 
of neutrinos, 

n ( ml \ / V2G F N e {x) \ 

H[6 cp ,N e (x)} = —\ m| C/t+ 0, (1) 

IK \ m 2 ) V / 

where E v is the energy of the neutrino, Gf the Fermi constant, and rrij, j=l, 2, 3 are the neutrino mass 
eigenvalues. U is the three-neutrino mixing matrix in the basis where the charged leptons are diagonalized Jl^ , 
|l5| and 5 cp is the CPV phase appearing in the mixing matrix. N e (x) is the DEN function that determines 
the matter effect. For the antineutrino N e (x) is replaced by —N e (x) and U by its complex conjugated which 
is equivalent to replacing S by —5. The v a — > vp oscillation probability and that for their anti-particles can 
be written succinctly to exhibit the DEN dependence as 

P a p = P Va ^ v JL,E, S cp , N e (x)) = \[Texp( - i I H[S cp , N e (x)]dx)) I, 

P a a = Pu a ^u {L,E,8 cp ,N e {x)) = \(rexp(-i( H[-6 cp ,-N e {x)]dx)) I, 

V J /a/3 
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where T denotes a propagation path-ordering product fljq . For numerical implementation we follow the 
usual approach of numerically integrating the Schrodinger equation. 

Although EDN is a critical factor in the analysis of the long baseline oscillation data, what is usually 
available is an averaged Earth density function N e (x) such as the widely used PREM model. In the following 
we study the effect of variations from the average value. Let us define the average density N e (x) and its 
uncertainties given by the variance function a(x), 

N e (x) =< N e (x) >= f[VN e }N e (x)F[N e (x)}, < a(x) >= yj '< N 2 e {x) > -< N e (x) > 2 , (3) 



where F[N e (x)][DN e ] is the probability of obtaining the DEN N e (x) in the neighborhood x, which will be 
defined in more detail later. The oscillation probability can be defined as an average over all the possible 
DEN profiles around N e (x). The appropriate framework for such a statistical expectation is the functional 
integral in which the Earth density profiles span a functional space that contains all possible variations of 
the earth densities allowed within the given variance. We write, 

< P a p >= J pN e }P af3 F[N e (x)}. (4) 

In geophysics, a uniform random sampling in a broad space of the Earth matter density are used to 
generate the matter density models ]l7, 18|. The density function generated in this way is further constrained 



by testing against two important sets of observational data. One is the mass and moment of inertia of the 
Earth, and the other the normal modes of the free oscillation of the Earth. In general the samples remained 
approach a Gaussian-like statistics rather than a uniform one, which together with the fact that the Earth 
matter density is always positive suggest to us to recast it into a logarithmic normal distribution for the 
probability density p9[, 

F[N e (x),x] = — -L=- exp(-\n 2 (N e (x)/N (x)) / (2s 2 (x))) 
N e (x)v2'Ks(x) v 7 

s(x) = yj\-a(l + r 2 {x)), N (x) = N e {x)exp( - s 2 {x)/2) (5) 

where r(x) = <j{x) / N e (x) parameterizes the uncertainty in DEN in terms of the ratio of local variance and 
the local mean value of the Earth's density. 

The logarithmic distribution is not a symmetric distribution for arbitrary a. However, it is close to the 
Gaussian distribution for a small in comparison with N e . In Fig. 1 we plot the logarithmic and Gaussian 
distributions for r = 5%, 50% respectively with N e (x) given by PREM. We see that the logarithmic distri- 
butions N e (x) is always positive and the difference between the Gaussian and the logarithmic distributions 
is very small in the case of r = 5% . In Fig. 2, we plot the density profile, from which we see that the PREM 
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denoted by the thick solid line is a average of the geophysical density samples shown by the oscillating thin 
lines. 

With the inclusion of uncertainties in DEN, we can now estimate the uncertainty in the neutrino oscil- 
lation probability by computing the variance of the oscillation probability 



8P a p = \]< {P afi - < P a/3 >) 2 > = y j {P a p- < P a! i >) 2 F[N e (x)}[DN e }. (6) 

The variance 5P a p leads to an uncertainty in the number of observed charged leptons. Since the charged 
lepton events are usually divided into energy bins, we define the variance of event number in a bin as, 

5N (E, L) = c(> Ua {E, L)5P a0 ap{E)AE, (7) 

where <p Va (E,L) is the neutrino beam flux spectrum of flavor a, crp(E) is the charged current cross section 
of neutrino of flavor /3, and AE is the bin size. 

To measure the CPV effect, the difference between event rates of opposite charged leptons is usually 
considered, 

N cp {E,L,S cp ) = Np-Np = [<t> Va (E,L)P afi a p (E)-<l> Pa (E,L)P a pa0(E)]AE 

= 4>,AE,L)D^(S cp )<jp(E)AE, (8) 

where the CP-odd difference |2(], ||, |J] is given by, 

D^(6 cp ) = P Va ^{L,E,8 cp ) - P Pa ^(L,E,S cp ). (9) 

In our numerical calculations, we have assumed that the neutrino and anti-neutrino beams have the same 
flux spectrum and the mass of the detector for the anti-neutrino is twice of that of the neutrino so as to 
compensate the difference in the neutrino and anti-neutrino charged-current cross sections. 
Now the error from the uncertainty of DEN is estimated as the standard deviation, 



5D c a ? = J < [P aP - P &f3 - < P afj - P a/ 3 >] 2 > 



[SP a p] 2 + [SP &0 ] 2 . (10) 

Unless this uncertainty is under control it will be difficult to extract the CP phase. It may even be difficult 
to establish the CPV effect if the fluctuation caused by the uncertainty of the matter density is not much 
smaller than the typical CP asymmetry given by, 

AD% = D aP (5 cp ) - D al3 (5 cp = 0). (II) 

Note that we can also work with the normalized conventional CP asymmetry, 



D CP / r 

ACP = ~5 — V~5 — ' 6ACP = \ / ( ACP ~ < ACP >fF[N e (x)][VN e ], (12) 

Pa/3 +P & f3 \ .1 

which should provide the same information. 

We evaluate Eqs. (4) and (6) numerically using a method similar to that of the lattice gauge theory. The 
neutrino path is discretized into / one-dimensional cells where / is a sufficiently large integer. In each of the 
?-th cell, the DEN function N e (xi) has an independent logarithm normal distribution with its local variance 
a(xi). Then the estimators of the mean and the deviation can be recasted respectively into the following 
forms, 



< 



D%> = Urn J D^[{Am^9;S cp } ] {N e {x 1 ),...N e {x i ),...N e (x I )}]l[F[N e (x i ) 

i—1 



dN e ( Xi ) 



N e (xi)V2Trs(x i ) 



= lim K-^bkiiAm 2 ,e,5 cp };{N e } k ] (13) 

fe=l 

K 1/2 

5D C J = { lim (K-l)- 1 J2iD k -<D^>} 2 } (14) 



K- 

k=l 
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where we have replaced the functional integration over the DEN by a sum over K arrays, {N e }k, k = 1, 2, ...K, 
and -Dfc is the value of the difference D^S over the /c-th density function {N e } k . The array, {N e } k which 
consists of density function N e (xi), ...N e {xi)...N e {xi) is generated from PREM weighted with a Gausian- 
like logarithm deviation. Of course, other specific Earth density models can be used. We have checked 
numerically that Eqs. ([il]) and (^) are convergent and stable. 

Geophysically, to obtain the density profile, the earth is discretized into hexahedrons ( elementary volume 
in spheroidal coordinates ), then Earth density is defined on the nodes and solved as an inverse problem. 
Generally there exists a length scale of the hexahedrons, which typically is order of 100 Km[jl0[ [HI ^J. So 
we take I to be the order of L/lOOKm in our calculation. Furthermore since the series in Eqs. (13) and (14) 
converges very rapidly due to the Gaussian-like distribution of Eq. (5), for the convenience of the calculation 
we identify K with the number of beam neutrinos in the individual bins, i.e., 

K = K(E, L) = 4> Va (E, L)a a (E)AE. (15) 

Having constructed all the needed ingredients, we can now investigate the error in CPV due the uncer- 
tainty of matter density. In the numerical calculation we take the baseline 2900 km for illustration. This 
baseline has been widely used for the study of CP violation at neutrino factories. We take a 20 GeV high 
performance neutrino factory, which delivers 10 21 working muons per year to a 50 kiloton detector. Since 
the path of the 2900 km baseline can go as deep as 160 km into the Earth and this will reach part of the low 
velocity zone, the uncertainty of the PREM can be sizable. In the numerical calculation, we take r(x) = 5% 
and the lattice size to be 200 km. 

We adopt the large mixing angle (LMA) scenario J2^] for the solar neutrino puzzle and take the following 
typical set of mixing parameters, 

Am 2 sol = 6.0 x 10~ 5 eV 2 , Am 2 tm = 3.55 x 10~ 3 eV 2 , 

tan 2 6 12 = 0.3, sm 2 26» 23 = 0.99 . (16) 

We also take the typical value sin 2 0i^ — 0.08. 

Now we present the numerical results. First, we study the extend to which the matter effect will mimic 
the CPV effect. In Fig. 3, we plot N cp as a function of the neutrino energy. The dashed line is given by 
S C p = with the error bars for an 5% uncertainty in PREM. So at this baseline the uncertainty caused 
by that of the matter density seems to be partially controllable. To see this more clearly and to estimate 
the range of the uncertainty, we show in Fig. 4 the following cases: dashed curve for S cp = 90° with 5% 
uncertainty in PREM, dotted and solid curves respectively for S cp = 54° and 5 cp — without uncertainty in 
PREM. It is clear that S cp — 90° can be easily distinguished from 5 cp = 0. However, there is an large error 
in extracting the CP phase that is about 36°. 

The uncertainty increases with the baseline because of the accumulation of the matter effect. In Fig. 5 we 
plot N cp vs the neutrino energy for a 12000 Km baseline. We consider again a neutrino factory which delivers 
10 21 working muons per year but at 50 GeV to increase the statistics. One can see that the uncertainty is 
large and it can no longer to distinguish S cp — 90° from S cp — 0° . So in order to have sensitivities for the 
CP measurement at this distance, the accuracy in DEN has to be much better than 5%. 

We note that the effect of uncertainties in the Earth density have also been examined in Ref. |24j. 
However the uncertainty has been fixed at the maximum value, i.e., N e (x) = (1 ± r')N e (x) where N e (x) 
is also given by PREM. As a comparison, we show in Fig. 5, together with our results, the constant 
uncertainties with r' = 5% (dot-dashed and dashed lines). In our language this fixed- value distribution 
function is represented by a delta function. So we think it overestimates the effects of the uncertainty of the 
Earth density. 

To conclude, we have considered in some details the issue of extract the information of the CPV effect 
in the presence of uncertainties in the Earth matter density in high precision measurements anticipated in 
future LBL neutrino experiments. We have developed a path integral formulation to estimate the fluctuation 
around the statistical mean. We have also presented a numerical implement of the formulation and applied 
it to assess the effectiveness in the determination of the CP phase. 

We thank Yi-Fang Wang for discussions. We also thank Fu-Tian Liu in the National Geological and 
Geophysics Institute for geophysics discussion. The work is supported in part by the NSF of China under 
Grant No 19925523 and also supported by the Ministry of Science and Technology of China under Grant No 
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Figure 1: Plot of probability distribution F[N e (x)][DN e (x)] vs density N e (x) at depth, e.g., R — 6200 Km 
from the earth surface. The mean value is taken from PREM which is Zg/cm? . The solid and dotted lines 
are the logarithmic and the Gaussian distribution with r = 5%, which are almostly coincided. The dash and 
dot-dashed lines are the logarithmic and Gaussian distribution for r = 50% respectively, which show that 
Gaussian distribution leads to negative density. 
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Figure 2: Plot of earth matter density vs path length. The thick solid line is the prediction of PREM 
N e (x). The dashed line above PREM is for N e (x) = (1 + 5%)N e (x). The oscillating lines are sample profiles 
generated according to Eq.(||) with r(x) = 5% . 
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Figure 3: Plot of CP-Odd event number difference N cp vs neutrino energy E v for baseline L=2900Km. The 
dashed line is for S cp = 0° and the error bars represents the variance caused by the uncertainty of matter 
density with r — 5%. The solid line is the prediction of S cp = 7.5° in PREM, i.e. r(x) = . 
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Figure 4: The same plot as Fig. 3, but parameters chosen differently. The dashed line is for 5 cp = 90° with 
error bars representing the uncertainties of r = 5%. The solid and dotted line are for S cp — 54° and 0° 
respectively with r = 0. 
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Figure 5: The plot of CP-Odd event number difference N cp vs the neutrino energy E v for baseline L — 12000 
Km. The dashed line is for S cp — 0° with error bars representing r — 5%. The solid line is for S cp — 90° with 
r = 0. The dash-dotted and dotted lines are given by N e (x) = (1 ± 0.05)./V e (a;) respectively, where N e (x) is 
given by the PREM. 
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